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A current objective of low-energy nuclear theory is to build non-empirical nuclear energy density 
functionals (EDFs) from underlying inter-nucleon interactions and many-body perturbation theory 
(MBPT). The density matrix expansion (DME) of Negele and Vautherin is a convenient method 
to map highly non-local Hartree-Fock expressions into the form of a quasi-local Skyrme functional 
with density-dependent couplings. In this work, we assess the accuracy of the DME at reproducing 
the non-local exchange (Fock) contribution to the energy. In contrast to the scalar part of the 
density matrix for which the original formulation of Negele and Vautherin is reasonably accurate, 
we demonstrate the necessity to reformulate the DME for the vector part of the density matrix, which 
is needed for an accurate description of spin-unsaturated nuclei. Phase-space averaging techniques 
are shown to yield a significant improvement for the vector part of the density matrix compared 
to the original formulation of Negele and Vautherin. The key to the improved accuracy is to 
take into account the anisotropy that characterizes the local-momentum distribution in the surface 
region of finite Fermi systems. Optimizing separately the DME for the central, tensor and spin- 
orbit contributions to the Fock energy, one reaches a few-percent accuracy over a representative set 
of semi-magic nuclei. With such an accuracy at hand, one can envision using the corresponding 
Skyrme-like energy functional as a microscopically-constrained starting point around which future 
phenomenological parameterizations can be built and refined. 
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I. INTRODUCTION 

The nuclear energy density functional (EDF) approach 
is the many-body method of choice to study medium- 
mass and heavy nuclei in a systematic manner [T] . Mod- 
ern parameterizations of empirical energy functionals 
(e.g. Skyrme, Gogny or their relativistic counterparts) 
provide a fair description of bulk properties and certain 
spectroscopic features of known nuclei. However, such 
empirical EDFs lack predictive power and a true spectro- 
scopic quality away from known data. Consequently, an 
intense ongoing effort is dedicated to empirically improv- 
ing the analytical form and the fitting of energy density 
functionals [2HU- 

A complementary approach in the quest for predictive 
EDFs [SHU] relies less on fitting empirical functionals to 
known data, but rather attempts to constrain the analyt- 
ical form of the functional and the values of its couplings 
from many-body perturbation theory (MBPT) and the 
underlying two- and three-nucleon (NN and NNN) inter- 
actions. Switching from conventional hard-core poten- 
tials to low-momentum interactions is essential in this 
respect, as the many-body problem formulated in terms 
of the latter becomes significantly more perturbative 1 . 



Indeed, second-order perturbative calculations provide a 
good account of bulk correlations in both infinite nuclear 
matter [T3] and doubly-magic nuclei [13]. Using many- 
body perturbation theory (MBPT) .15] as a baseline, 
the long term goals of the project are to (i) bridge non- 
empirical EDF methods with ab-initio many-body tech- 
niques applicable to light nuclei, (ii) calculate properties 
of heavy/complex nuclei from basic vacuum interactions 
and (iii) perform EDF calculations with controllable the- 
oretical errors. 

MBPT contributions to the energy are written in 
terms of density matrices and propagators convolved with 
finite-range interaction vertices, and are therefore highly 
non-local in both space and time. In order to make such 
functionals numerically tractable in heavy open-shell nu- 
clei, it is desirable to develop simplified approximations 
expressed in terms of the local densities and currents. 
Starting at lowest order, which displays only non-locality 
in space through the Fock contribution to the energy 2 , 
the objective of the present work is to revisit the density 
matrix expansion (DME) of Negele and Vautherin jTHj to 
assess its accuracy in reproducing non-local Fock contri- 
butions. 

The focus of the present paper is on the vector part of 
the density matrix, which is relevant for approximating 
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1 The need for infinite resummation of certain sets of diagrams 



and/or the redefinition of the unperturbed vacuum |<J>) cannot 
be ruled out at this point. 
2 For simplicity, we are assuming local NN and NNN interactions. 
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EDF Energy density functional 
DME Density matrix expansion 
PSA Phase space averaging 
OBDM One-body density matrix 
INM Infinite nuclear matter 



TABLE I: List of acronyms repeatedly used in the text. 



perturbed many-body state around which perturbation 
theory is performed or, in a more phenomenological 
language, the auxiliary state in terms of which one 
builds a so-called single-reference energy density func- 
tional (EDF). In the present case, we consider an im- 
plementation without explicit treatment of superfluidity 
such that |$) takes the form of a Slater determinant. In 
addition, we consider the system to be invariant under 
time-reversal. 



the central, tensor and spin-orbit Fock contributions in 
spin-unsaturated nuclei, i.e. in nuclei where only one of 
two spin-orbit partners is filled. Indeed, the few tests 
of the DME over the past thirty-five years have focused 
entirely on the scalar part [TU [TTl [18] , given that no re- 
liable expansion of the vector part of the density matrix 
was ever proposed. As acknowledged by Negele and Vau- 
therin in their seminal paper, the expansion suggested 
for the vector part of the density matrix was not on the 
same level as the one designed for its scalar part. Such a 
feature is obviously critical since the overwhelming ma- 
jority of nuclei are spin unsaturated. Here, we demon- 
strate that phase-space averaging techniques allow a con- 
sistent expansion of both the scalar and the vector parts 
of the density matrix, such that the accuracy is greatly 
improved for the latter. A key feature of the new method 
is to take into account the deformation displayed by the 
local momentum distribution at the surface of most fi- 
nite fermi systems [THl [2H]. While it is shown to have 
little impact on the expansion of the scalar part, the de- 
formation of the local momentum distribution is crucial 
to accurately reproduce contributions to the energy that 
probe the vector part of the density matrix. 

The paper is organized as follows. Section |H| provides 
the basic ingredients needed to conduct the present study. 
Section Mil is dedicated to the reformulation of the den- 
sity matrix expansion on the basis of phase-space aver- 
aging techniques. The accuracy of the approximation 
method is gauged in Sec. |IV| through non-self consistent 
tests that make use of two schematic nucleon-nucleon 
interactions and of density matrices obtained from self- 
consistent EDF calculations of a large set of semi-magic 
nuclei. Each of the central, tensor and spin-orbit contri- 
butions to the Fock energy is analyzed separately. Con- 
clusions are given in Sec. [V] while appendices provide 
complete sets of formulae and analytical derivations. In 
particular, couplings of the generalized Skyrme-like EDF 
obtained through the DME (see Eq. 20 ) are provided in 
appendix [5] 



II. DENSITY MATRIX AND HF ENERGY 



A. The one-body density matrix 

The one-body density matrix (OBDM) p of the 
many-body state |<fr) is defined in terms of operators 

(r a q)/c (r a q) that create/annihilate a nucleon at a 
given position in space r with given spin and isospin pro- 
jections cr = ±1/2 and q = n,p on the quantization axis 

p q (ra,r'a') = (<f>| c f (r' a' q) c (r a q) |$) 



where it is assumed that single-particle states do not 
mix isospin projections so that the OBDM is diagonal in 
isospin space 3 . In Eq. [TJ p^ = ( $ | cj Cj | $) defines the 
OBDM in an alternate single-particle basis {c,;; ipi(r<rq)}. 
Choosing the particular basis from which |<E>) is built, 
becomes diagonal with matrix elements equal to one for 
occupied states and zero for empty states. The OBDM 
can be further separated into 



p q (ra,r'a') = ^{p g (r,r')<W + s g (r,r') • cr CTCT '} ,(2) 



where the scalar and vector parts are respectively defined 
as 



p q (r,r') = ^p 9 (r<7,rV)<</|l|a) 

era' 

= mi ^(r'^^rcrg) 

a ij 



(3) 



= EE^^V^^Vl^^-^p], .(4) 

aa' ij 

In the approximation that the single-particle wave- 
functions of spin-orbit partners are identical, it can be 
shown that the vector part of the density matrix s g (r, r') 
is zero in spin-saturated nuclei. 



Let us consider a product state of reference |$). As 
briefly explained in Sec. II B| this typically is the un- 



3 The Slater determinant can however break spatial symmetries. 
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B. Long-term strategy 

Our long-term objective is to build so-called non- 
empirical nuclear energy functionals E [p] through the ap- 
plication of many-body perturbation theory implemented 
in terms of low- momentum interactions 121 



£[p] =E Hh +AE 



(5) 



where E HF denotes the (symmetry-unrestricted) 
Hartree-Fock (HF) contribution from two-, three- 
. . . nucleon forces whereas AE HF encompasses the 
corresponding correlation energy to all orders in per- 
turbation theory 4 . As opposed to the wisdom based on 
the use of conventional nuclear potentials, it has been 
shown recently that so-called low-momentum two- and 
three-nucleon interactions make the nuclear many-body 
problem more perturbative, with Hartree-Fock serving 
as a reasonable zeroth-order approximation |13j . Still, 
calculations of the infinite nuclear matter (INM) equa- 
tion of state [15] . as well as binding energies and charge 
radii of doubly- magic nuclei [14j . demonstrate that it is 
necessary to go at least to second-order in perturbation 
theory to resum enough bulk correlations into the EDF 
to get realistic binding. In the present paper though, 
we focus on the lowest-order contribution to the energy 
that is bilinear in the OBDM, i.e. the Hartree and Fock 
diagrams. While treating the direct (Hartree) term 
exactly, the objective of the density matrix expansion 
is to simplify the non-local character of the exchange 
(Fock) contribution to the energy by mapping it into a 
generalized Skyrme functional with density-dependent 
couplings. Therefore, the DME can be viewed as a 
constructive approach to encode finite-range physics into 
density-dependent couplings of a Skyrme-like functional. 

The reasons for restricting our attention to the 
Hartree-Fock contributions in this initial study are two- 
fold. First, a non-trivial extension of the DME is needed 
to treat non-localities in both space and time that arise 
in higher orders of perturbation theory. I.e., one must 
properly account for the presence of energy denominators 
when designing a DME for 2nd-order MBPT and beyond. 
To date, a satisfactory generalization of the DME has not 
yet been formulated. Second, even if we follow the ad- hoc 
prescription of neglecting the non-locality in time by us- 
ing averaged energy denominators, it is well established 
that the dominant contributions to bulk nuclear proper- 
ties are of the Brueckner-Hartree-Fock (BHF) type. Op- 
erationally, this amounts to replacing the vacuum NN in- 
teraction in the Hartree-Fock expression by a Brueckner 



4 In applications to nuclei, except for doubly-magic ones, the 
ground-state energy will in fact be expanded around a quasi- 
particle vacuum of the Bogoliubov type rather than around a 
Slater determinant. This is necessary to take care of the Cooper 
pair instability that arises in the ^-So channel of the in-medium 
NN amplitude. 



G-matrix (or a perturbative approximation in the case 
of low-momentum interactions) evaluated at some aver- 
age energy. Since the G-matrix "heals" to the NN po- 
tential at long distances, applying the DME to the long- 
range part of the NN interaction at the Hartree-Fock level 
will in any event capture the same contributions to the 
density-dependent couplings as given by the long-range 
part of the G-matrix in a more sophisticated BHF cal- 
culation. In this way, the dominant density-dependence 
that arises from the finite-range of the inter-nucleon in- 
teractions is accounted for. Once a satisfactory general- 
ization of the DME is developed to handle spatial and 
temporal non-locality on the same footing, non-localities 
arising from in-medium propagation can be mapped into 
the density-dependent Skyrme couplings as well. 



C. Two- nucleon interaction 

For simplicity, and because the main point of the 
present paper does not depend on it, we restrict our study 
to two-nucleon interactions only. Note however that a 
forthcoming publication is dedicated to the application of 
the presently developed DME to the HF energy derived 
from a chiral-EFT three-nucleon potential at N 2 LO |22) . 
In the present paper, we consider a generic local two- 
body interaction that includes central, tensor and spin- 
orbit parts. Defining Xi = (r^cr^), one can write in the 
position ® spin ® isospin basis 

(xix 2 \Vf T \x 3 x 4 ) = Vf T 5(n - r 3 ) <5(r 2 - r 4 ) ,(6) 

where / can be G— central, LS— spin-orbit or T— tensor 
whereas (S,T) takes values (1, 0), (0, 1), (1, 1) or (0,0), 
where the first number 1/0 refers to two-body spin- 
triplet/singlet channels whereas the second number 1/0 
refers to two-body isospin-triplet/singlet channels. More 
explicitly, the central part of the interaction reads 



V, 



ST 
C 



ST 

J C 



(r) n; 



It L1 s/t 



where the relative and center of mass coordinates are 
defined as 



r = ri - r 2 and R = ^( r i + r 2)- 
while spin/isospin singlet /triplet projectors 



H 



s/t 



-(1 -/+ Ff 2 ) and W sJt = -(1 -/+ PI, 



(7) 



(8) 



are expressed in terms of spin/isospin exchange operators 



Pl2 = \{^2 



1) and P[ 2 = l( n .T2 + l). (9) 



The spin-orbit and tensor parts of the two-nucleon inter- 
action take the form 



'LS 



ST _ „,ST 



(r) 



3(cti • e r ) (cr 2 • e r ) — o\ ■ cr 2 
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with e r = r/r. It should be noted that the spin-orbit and 
tensor parts of the interaction only act in the spin-triplet 
channel. 



D. Fock contribution to the energy 

As mentioned earlier, the strategy consists of applying 
the DME to the exchange part of the HF energy while 
treating the Hartree term exactly. Indeed, it was realized 
long ago, starting with the early works on the DME by 
Negele and Vautherin [TH] H3J, that treating the direct 
part exactly has the following advantages: 

(i) It provides a better reproduction of the density fluc- 
tuations and the energy produced from an exact HF 
calculation [25] . 

(ii) It significantly reduces the self-consistent propaga- 
tion of errors if one restricts the DME to the ex- 
change contribution [TTJ [23] . 

(hi) There is no additional complexity in the numerical 
solutions of the resulting self-consistent HF equa- 
tions [23] compared to applying the DME to both 
Hartree and Fock terms. 

The Fock contributions from central, spin-orbit and 
tensor parts of the two-body interaction take the form 



E£[ST] 



dr\dr 2 



/3 g (ri,r 2 )p g '(r 2 ,ri) 



s g (ri,r 2 ) • s g /(r 2 ,ri) v^ T (r) , (10) 



drirfr 2 



p g (ri,r 2 )r x V 2 • s ? /(r 2 ,ri) 



s g (ri,r 2 ) -r x V 2 pg'(r 2 ,ri) v[^(r), (11) 



such that the exchange contribution from the spin-orbit 
interaction reduces to 



E[ S [ST] ~ J dr 1 dr 2 vf§(r)s q {r 1 ,r 2 ) -r x V 2 p q >(r 2 ,r x ) 
III. REVISITING THE DME 



Basics of the DME 



The DME was originally proposed by Negele and Vau- 
therin to establish a theoretical connection between the 
empirical zero-range Skyrme force and Hartree-Fock cal- 
culations with realistic NN interactions [T5] . The central 
idea is to factorize the non-locality of the OBDM by ex- 
panding it into a finite sum of terms that are separable 
in relative and center of mass coordinates. Adopting no- 
tations similar to those introduced in Refs. [23] and [26] . 
one writes 



P ,(n,r 2 ) » Y^nr n (kr)V n {K) , (15) 

m max 

s„(ri,r 2 ) « J2 H ™( fcr ) 2 ™( R ) > ( 16 ) 

m=0 



where k is a momentum scale to be determined that 
sets the scale for the decay in the off-diagonal direc- 
tion, n^(fcr) are the so-called n~ functions that re- 
main to be specified, and {"P„(R), Q m (R)} denote var- 
ious bilinear products of local densities and their gradi- 
ents {/j g (R),r 9 (R), Jq, MV (H), Vp g (R), Ap q (R)} obtained 
from the OBDM through 



E$[ST\ 



dridr 2 



s g (ri,r 2 ) • s ? /(r 2 ,r!) 



+ E^ s ?>i< r 2K>(r 2 ,ri)] v s T T (r)\12) 



where numerical coefficients and overall signs, as well as 
sums and/or selection rules over isospin projections have 
been omitted. Indeed, only the structure of the terms at 
play is of importance for the present paper. For time- 
reversal invariant systems, the scalar and vector parts of 
the OBDM satisfy the relations [23] 



p q (r x ,Y 2 ) = p q {Y 2 ,r 1 ) , 
s 9 (ri,r 2 ) = -s 9 (r 2 ,ri) 



(13) 
(14) 



p q (R) = p 9 (ri,r 2 )| ri=r2= R , 
T q (R) = Vi • V 2 p g (r 1 ,r 2 )| ri , 



= R 



(17) 
(18) 



~^1~^2)„ M r l> r 2)L=r 2 =R (19) 



The above local densities relate to the matter density, 
the kinetic density and the cartesian spin-current pseu- 
dotensor density, respectively. See Appendix [X] for more 
details. Provided that large enough n max and m max give 
an accurate reproduction of the Fock contributions to the 



energy (Eqs.llO 11 and 12 1, the benefit of expansion 15 



16 is to provide a local approximation of the form (for 
time-reversal invariant systems) 
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E F ^J2 I dnlAPP Pq {H) Pq {H) + A? T p q (R) Tq (R) + A? A » p q (R) Ap,(R) + A<> VJ p q (R)V • J,(R) 

q J < 

+ A^P J W Pq (R) ■ 3 q (R) + A JJ Y, -W R ) J q ,^CR) 

fjLV 

+ A JJ [fe J 9.mm(R)) (E-WR)) + E } 
+ E / rfR{s pp p 9 (R)p,-(R) + B" r p 9 (R)r 9 -(R) + Br A »p q (R)A Pq (R) + B pvJ p 9 (R)V ■ J 9 -(R) 
+ B v P J W Pq (R) • J ff (R) + 5 JJ J] J fl)/ ,„(R) Jg^(R) 

E ) ( E ^3.w( R ) ) + E ^:Aw( R ) ^^(r) 



(20) 



which is nothing but a local Skyrme-like EDF with cou- 
plings microscopically derived from the vacuum interac- 
tion. The couplings depend on the yet to-be-specified 
momentum scale k, and are given by integrals of the 
finite-range NN interaction over various combinations of 
Il-functions, e.g. 



A pp [fc]-47r r 2 drv^(r) ng(fcr) 



(21) 



Complete formulas for all the couplings appearing in 
Eq. [20] are provided in appendix [Bj Before coming to the 
details of the DME method, a few remarks are in order: 



(i) Eventually, the momentum scale k will be linked 
to the local Fermi momentum k F (R) , or to a simi- 
lar function, such that all couplings become den- 
sity/position dependent. 
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one sees 



From Eq. 

that such density /position dependence is a direct 
consequence of the finite-range of the NN interac- 
tion 



20 



is 



In this respect, the form given in Eq 
more general than any existing empirical Skyrme 
EDF. 

(ii) Due to such a density/position dependence of 
the couplings, terms that are usually connected 
through a partial integration, e.g. p q (R) A p q (R) 
and Vp g (R) • Vp g (R), can in general no longer be 
transformed into one another. As a result, one 
keeps both types of terms explicitly in the resulting 
EDF. 

(iii) Starting from a realistic vacuum Hamiltonian con- 
taining a three-nucleon force, one obtains a richer 
EDF including a wealth of trilinear terms |22| . In- 
cluding such terms will be eventually essential to 
any realistic application of the present work. 



iv) Eq. 20 is to be complemented with the Hartree con- 



tribution that can either be put under the form of a 



local EDF or treated exactly. Regardless, the EDF 
thus obtained only contains the physics of the HF 
approximation such that further correlations must 
be added in order to produce any reasonable de- 
scription of nuclei. In the short term, such an ad- 
dition can be done empirically by adding the above 
DME coupling functions to empirical Skyrme func- 
tionals and performing a refit of the Skyrme con- 
stants to data. This phenomenological procedure is 
motivated by the earlier observation that a Brueck- 
ner G-matrix differs from the vacuum NN interac- 
tion only at short distances. Therefore, one can in- 
terpret the refit to data as approximating the short- 
distance part of the G-matrix with a zero-range ex- 
pansion thru second order in gradients. Eventually 
though, and as already stated, it is the goal of a 
future work to design a generalized DME that is 
suited to higher orders in perturbation theory. 



B. Existing variants of the DME 

Several DME variants applicable to the HF energy 
have been developed in the past .16, 27 -29]. They mainly 
differ regarding (i) the choice made to fix the momentum 
scale k, (ii) the path followed to obtain actual expres- 
sions of the il-functions (see below) and (iii) the set 
of local densities that occur in the expansion. For in- 
stance, the DME of Ref. [37] is a variant of the original 
one proposed by Negele and Vautherin (NV-DME) [TB] 
that improves the accuracy of the expansion obtained at 
first order (n max = 0) by optimizing the momentum scale 
k. The DME of Ref. is based on a semi-classical ex- 
tended Thomas-Fermi approximation, while the one pro- 
posed in Ref. [28] is a phenomenological method that 
introduces parameters to be optimized in order to obtain 
the correct local semiclassical kinetic energy density and 
integrated projector identity of the OBDM (see Eq. 30 1. 
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C. Motivation for a PSA reformulation of the DME 

The central part of the present work relates to a new 
and more general DME variant that is based on phase- 
space averaging (PSA) techniques. It will be denoted as 
PSA-DME throughout. The need for such a new for- 
mulation of the DME, in light of the number of already 
available variants, relies on the following observations 

(i) Existing DME formulations have focused mostly on 
the scalar part of the OBDM. For instance, Negele 
and Vautherin acknowledge in their seminal paper 
that they were not able to design an approximation 
of the vector part of the OBDM on the same level, 
and thus with the same accuracy, as the one they 
obtained for the scalar part. This is an essential 
problem in view of constraining non-empirically the 
nuclear EDF. Indeed, the vector part of the OBDM 
is non zero in spin-unsaturated nuclei, i.e. in almost 
all nuclei. 

(ii) The PSA reformulation proposed below provides a 
consistent derivation of the DME expansion of both 
the scalar and the vector pieces of the OBDM. In 
addition, it recovers the NV-DME as a particular 
case, such that one is offered the freedom to choose 
in a consistent fashion the variant that best opti- 
mizes the reproduction of the each of the three Fock 
contributions to the energy. 

(iii) In the PSA approach, one uses information from 
the local momentum phase space distribution of the 
system of interest in order to optimize the DME 
length-scale k and to produce analytical expres- 
sions for the Il^(fcr) functions. 

(iv) Finally, it should be pointed out that all available 
DME techniques hold only for time-reversal invari- 
ant systems. Hence, an approach that can be ex- 
tended to non time-reversal invariant systems is 
important to constrain the nuclear EDF for non- 
time reversal invariant systems. In that respect, 
the requirements of Galilean, alternatively gauge 
invariance, can be used to establish various rela- 
tions between the n— functions multiplying certain 
time-even and time-odd densities [351 13S] ■ 

Note that the PSA formulation of the DME is not com- 
pletely new. Negele and Vautherin mentioned the possi- 
bility to use such an approach, having in mind to use the 
phase space of infinite nuclear matter, before reverting to 
a formal Bessel-function plane-wave expansion. From a 
formal point of view, the PSA approach developed below 
differs from that mentioned in Ref. [16 and is applied 
consistently to both the scalar and the vector parts of 
the OBDM. For instance, in spite of the weak angular- 
dependence of the scalar part of the OBDM [31], the 
inconsistency in the order of application of the angle- 
averaging and series expansion that exists in Ref. [16 is 
not an issue in the present case. 
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FIG. 1: (Color online) The quadrupole anisotropy P^ l (R) of 
the local neutron momentum distribution in a selected set of 
semi-magic nuclei. The black, red and blue vertical lines in- 
dicate the approximate half-radii (where the density becomes 
half of the density at the origin) . 



D. Momentum phase-space of finite Fermi systems 

A finite fermi system exhibits peculiar properties for 
the momentum phase-space distribution that are not 
present for homogeneous systems. The intent of this 
section is to mention those features that are relevant to 
the present work. The local momentum distribution of 
quantum systems can be studied via a multitude of quan- 
tum phase-space distribution functions |32) . Using the 
Wigner distribution in Ref. [19] and the Husimi distri- 
bution in Ref. |20j . the local single-particle momentum 
distribution is shown to display a diffuse and anisotropic 
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Fermi surface when sitting at the (spatial) surface of the 
finite system. For reasons discussed in Sec. |III E the dif- 
fuseness is not as important as the anisotropy. Hence, we 
now describe a method that can be used to quantify of 



the anisotropy of the local Fermi surface. 

In Ref. [20], the local quadrupolar deformation of the 
momentum Fermi surface (for a given isospin) is given 
by 5 



i?(r) 



fdp[3(e r -p) 2 -p 2 ]H q (r,p) 
fdpp 2 H q (r,p) 



r g (r) 

r 



^|(e r .Vfe(rg)| 2 ^-l 



+ 0((k F r a ) 2 ), 



(22) 



where H q (r,p) is the Husimi distribution, tq is a length 
scale used in the Husimi distribution and k q F is a short- 
hand notation for the local Fermi momentum k F (R) de- 
fined in a local density approximation through 



k F — 



3tt 2 Pq {K) 



1/3 



(23) 



Equation 22 is computed in the basis <pi(rq) that diago- 
nalizes p, i.e. the basis from which the Slater determinant 
|$) is built 6 . A simplified expression of P$ (r) in spher- 
ical symmetry suitable for semi-magic nuclei is provided 
in appendix [C] 

Fig. [T] shows the quadrupole anisotropy of the local 
neutron momentum distribution calculated for a selec- 
tion of semi-magic nuclei. Single-particle wave-functions 
are obtained from a Skyrme-EDF calculation performed 
with the BSLHFB code [32] using the SLy4 parametriza- 
tion of the Skyrme EDF with no pairing. Figure T] also 
displays the local neutron Fermi momentum (Eq 
order to locate the position of the nuclear surface 
spite of pronounced shell fluctuations, the result corrob 



23|) in 
In 



orates the conclusions drawn in Ref. [20] ; P£ (R) becomes 
negative just inside the surface, denoting an oblate mo- 
mentum Fermi surface while, outside this region, the lo- 
cal momentum Fermi surface becomes strongly prolate. 
In both cases, we have taken an axis normal to the nu- 
clear surface as the reference axis. The next two sections 
show how we make use of these properties of the phase- 
space distribution of finite Fermi systems to design our 
PSA-DME. 



E. The scalar part of the OBDM 

In a nutshell, the PSA approach consists of three basic 
steps: (i) the isolation of the non-locality as an expo- 
nential derivative operator acting on the OBDM, (ii) the 
expansion of that operator around a momentum scale k 
and (iii) the averaging of that momentum scale over the 
local momentum distribution of the system of interest. 

Applying the first two steps to the scalar part of the 
OBDM of a time-reversal invariant system, one writes 



/ r _ r 
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(24) 



Before approximating the action of the non-locality op- 
erator, e r '(Vi-v 2 )/2^ a p nase factor e lr k was extracted 
in order to perform a Taylor series expansion of the non- 
locality about the momentum scale k. We presently trun- 
cate the expansion at second order although nothing pre- 
vents to study higher orders in principle. The next step 
consists in performing an angle averaging over the orien- 
tation of r, which is a reasonable step as the scalar part of 



the OBDM has negligible dependence on the orientation 
of r [3T]. See appendix [D] for details. 

The final step involves averaging the dependence on 
the momentum scale k over a model phase space that 
characterizes the system under study. Performing the 
PSA of a function <?(k) over the locally-equivalent pure 
isospin infinite matter phase-space, i.e. defining G(k q F ) 



as one obtains for time-reversal invariant systems 

G(k%) = * / dkg(k) (25) 
47rfcJ, J\k\<k% 



ng(fc«r)p,(R) 



¥ ng(4r) 
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with 



ng(fc|,r) 
H^r) 



3i(k q F (R)r) 
k q F (R)r 

ji(4(R)r) 
*|.(R)r 



(27) 
(28) 



For details of the derivation, refer to appendix[Dj Several 
comments are in order: 

(i) The phase space of finite nuclei has a marked dif- 
ference from that of INM PH ED]- Still, using 
INM phase space suffices for the scalar part as 
will be apparent from the results discussed in sec- 



tion IV B This is because, unlike the vector part 



of the OBDM discussed below, the scalar part is a 
bulk quantity with most of its contribution coming 
from the interior of the nucleus where, to a good 
approximation, the momentum distribution resem- 
bles the one of INM [31]. 

(ii) Dealing separately with the neutron or proton 
OBDM in a finite nucleus, it is natural to perform 
the corresponding PSA over the phase space of the 
locally-equivalent neutron or proton infinite mat- 
ter. However, this provides II— functions with an 
explicit isospin dependence that eventually breaks 
the explicit isospin invariance of the EDF (but not 
its isospin symmetry). Considering the small dif- 
ference between k q F and the total local momen- 
tum kp(R), defined in terms of the total density 
p(R) = p„(R) + p p (R) through 



3tt 2 



p(R) 



1/3 



(29) 



it might be preferred to perform the PSA over the 
phase space of symmetric nuclear matter, even in 
a neutron rich nucleus. In any case, all results pre- 
sented below are obtained using k q F but would not 
be significantly different if using k F instead. 

(iii) The DME is not a naive Taylor expansion of the 
OBDM with respect to the non- locality r. The 
II— functions resum dependencies on r to all or- 
ders such that the long distance limit behavior of 
the OBDM is reproduced (see below). However, as 



noted in Ref. |16j , the truncation of the expansion 
about k to second order leaves the specific value 
of the coefficients of terms beyond k q F r undeter- 
mined (in the Taylor series expansion of n^/c^r)). 
This indeterminateness gives one the freedom to 
optimize n^, which can be viewed as selecting a 
different rearrangement and truncation of the ex- 
pansion [IS] . 

(iv) The zeroth-order n— function I\^{k q F r) found above 
is exactly the one found in the original NV-DME 
of Ref. [H>]. Just as in the NV-DME, the lead- 
ing term of the PSA-DME reproduces the exact 
OBDM of infinite nuclear matter. The second or- 
der n— function n^fc^r) is different 7 from the one 
found in Ref. [IB] . However, this relates to the 
previous remark that emphasized the freedom in 
choosing the second-order n-function. Moreover, 
we will find in Section Hvl that these differences are 
rather small for contributions to the Fock energy. 
Therefore, our PSA-DME of the scalar part of the 
OBDM is essentially equivalent to the NV-DME of 
Ref. pi]. 

The freedom mentioned above can be used to adjust 
to satisfy certain properties of the exact OBDM, or sim- 
ply to optimize the quality of the approximation through 
a comparison with realistic a OBDM. One example re- 
lates to the integrated idempotency of the OBDM, e.g. 
for neutrons 



N = Jdrp n (r) = J J dv x dv 2 |p„(ri, r 2 )| : 



(30) 



As shown in Ref. [S3] , there is a class of DME that sat- 
isfies this constraint. Unfortunately, the given in Eq. 
(28) does not satisfy this constraint. Even though the 
non-self consistent result given in |IV B| is satisfactory, 
this might not be the case in a self-consistent test. 

Other constraints on the n— functions come from the 
expected limits for large and small values of r. The 
n— functions should go to zero in the large r limit, while 



7 The Bessel expansion of Ref. |16| provides IIj 
105i3(fcf,r)/(fcf,r) 3 . 
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for small r, the expansion must to reduce to a simple 
Taylor series. These requirements 8 lead to [25l [30] 



ng(o) 



(0) = 



ttP' 

lim = 



n§(o) = i 

lim II2 = 



(0) 



(31) 
(32) 
(33) 



It can easily be shown that the above constraints are 



satisfied by the II— functions listed in Eqs. (27 1 and (28 1 



F. The vector part of the OBDM 

Restricting again the discussion to time-reversal invari- 
ant systems and applying the same steps as for the scalar 
part of the OBDM, one obtains for its vector part 
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(34) 



n=r 2 =R 



where only the first order term in the expansion of 
the non-locality operator was kept for reasons explained 
below. One also notes that the zero-order term pro- 
vides the local spin density s g (R) which is zero for the 
time-reversal invariant systems we are considering. In 
Ref. [IB] , it was argued that averaging over the orien- 
tation of k and setting k = k F should be sufficient to 
provide a reasonable account of the vector part of the 
exact OBDM. This gives 



R 



,R 



iUl(k q F r) 



where 



n?(4r) = JO (4(R)r) 



lJq,/J.v(R-) 1 (35) 



(36) 



If instead one applies the same procedure as for the 
scalar part of the OBDM and performs the PSA over 
the locally-equivalent pure-isospin infinite matter phase- 
space, one obtains 9 



Ul(k F r) 



, Ji(fc|-(R)r) 
k q F (R)r 



(37) 



However, as mentioned in section |III D the local mo- 
mentum distribution in the surface region of a finite nu- 
cleus has a markedly different behavior than the isotropic 
momentum distribution of infinite nuclear matter. Given 
that the vector part of the density matrix peaks around 
the nuclear surface, it seems more appropriate to perform 
the PSA over a deformed Fermi sea that incorporates the 



See appendix [e| 



for details. 



information contained in the function P q (R) discussed in 
section |III D[ The details are given in appendix [E] The 
final result differs from that in Ref. [IB] only in the ana- 
lytical form of Ilf . The result reads 



ni(k q F r) 



ji(*gjgOr) 
~k q F (H)r 



where 



k F — 



2 + 2P|(R) 

2-P|(R) 



1/3 



k q F (K) 



(38) 



(39) 



The PSA over the locally-equivalent neutron or proton 
infinite matter modifies the analytical form of ilf com- 
pared to NV-DME, i.e. compare Eqs. [37]and[38| In ad- 
dition, and contrary to the scalar part of the OBDM for 
which it is unimportant, taking into account the defor- 
mation of the local momentum distribution of the finite 
system leads to a modification of the relevant momentum 
scale k F . In view of isolating the significance of such an 
effect, while preservin g the benefit of usin g PSA, one can 
set P|(R) 



in Eq. 39 In Sec. IV C 



we discuss and 
compare the accuracy obtained using all of the preceding 
variants. 

Note that the expansion was limited to first order in 
Eq. |34| The reason is that, for time-reversal invariant 
systems, the cartesian spin-current pseudotensor density 
Jg iAI „(R) and its gradients are the only standard local 
densities at hand to express the DME. Given that, we 
could not find any closed and parameter-free expression 
of higher-order contributions in terms of such local den- 
sities only. This points however to the possibility to 
study higher-order terms in the context of the generalized 
Skyrme EDF discussed in Ref. [B]. 
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Finally, one can easily verify that the large and small 
r limits, viz, 



n^(0) = 1 , 115 (0) = and lim Il s = 



(40) 



mentioned at the end of section |III E| are satisfied by the 
expressions of Ilf given by cither Eq. (37) or Eq.(38). 



IV. COMPARING PSA- AND NV-DME 

The accuracy of our newly developed PSA-DME needs 
to be tested against both non-self consistent and self- 
consistent HF calculations. A self-consistent test of the 
PSA-DME is the aim of a forthcoming publication. As 
explained below, we limit ourselves in the present paper 
to gauging the accuracy of the NV-DME and the PSA- 
DME against two non self-consistent measures. Where 
relevant, we also set Pf (r) = in the PSA-DME of the 
vector part of the OBDM to isolate the significance of 
using a deformed local momentum Fermi surface. We 
denote that last variant as INM-DME. 



A. Inputs to non-self-consistent tests 



independently of the (S,T) channel and with vq = 50 
MeV, a = 1.5 fm, m n = 0.7 frn . The momentum cut- 
off A is set equal to 2.1 fm -1 while erfc is the comple- 
mentary error function. It must be stressed that none of 
these interactions are realistic two-nucleon interactions, 
but rather schematic representatives. The objective of 
the present study is to gauge the accuracy of the DME 
variants against a reasonable reference point that is not 
itself meant to provide useful or realistic results. The 
application of the present DME scheme to realistic chiral 
two- and three-nucleon interactions is the objective of a 
forthcoming publication 22 . Finally, note that neutron 
density matrices and local densities used in the following 
sections have been obtained, for all semi-magic nuclei of 
interest, through spherical self-consistent EDF calcula- 
tions employing the SLy4 EDF parameterizations with 
no pairing. 

B. Fock contribution from Vc 



The expression of the Fock contribution to the energy 
from the central part of the two-nucleon interaction is 
given in Eq.(10). It contains a bilinear product of non- 



local matter densities as well as a bilinear product of 
non-local spin densities. Since the latter also appears as 
part of the tensor contribution to the Fock energy (see 
Eq.(12)), we postpone the discussion regarding the spin- 



The generic form of the central, spin-orbit and t ensor density product to section [IVC 



interactions considered here have been given in Sec. II C 



The radial form factors used in the present calculations 
for either of those interactions take the form (i) a gaussian 
or (ii) a renormalized Yukawa (according to Ref. [55]). 
Specifically we use 



v e 



Before comparing the Fock energy to its DME counter- 
part, we first conduct a more stringent test on the energy 
density in which the integration over the angle of r has 
already been performed, i.e. we compare the integrand 



Cn«( R > r ) = ~ Jde rPn (r 1 ,r 2 )p n (r 27 r 1 ) , (42) 



[ e - m - r erfc(^ - rX) - (r -» -r)] , 



(41) to its DME counterpart 



a 



DME 



Rr) = ITg(fc£r) p n (R)p n (R) + ~ng(fe£r)n£(fc£r)/>„(R) 7 Ap„(R)-r„(R) 



jjfc£ 2 p„(R) 



,(43) 



where the latter depends on which variant of the DME 
has been adopted 10 . Having in mind existing empirical 
Skyrme EDFs that contain only up to two spatial deriva- 
tives, terms containing fourth-order gradients have been 



10 We denote such integrands as energy densities throughout the 
paper. Strictly speaking, it is necessary to multiply them by the 
interaction to obtain the dimension of an energy density. Still, we 
postpone the folding with the interaction to the second measure 
introduced below. 



truncated in C^' 1E (Il, r). A consistent account of such 
fourth-order derivatives in the EDF would require to go 
also to fourth order in the DME itself, which is beyond 
the scope of the present study. This is an important point 
that underlines our philosophy that the primary purpose 
of the DME method is not to reproduce the fine details of 
the OBDM, but rather to reproduce as best as possible 
the energy density and the total energy at a given order 
in the expansion. The latter two are precisely what is 
gauged in this paper, whereas no tests dedicated to the 
reproduction of the OBDM by itself are performed. 
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do here. Note also that, although the plots are provided 
for two sample nuclei, more systematic tests have been 
performed over several semi-magic isotonic and isotopic 
chains that support such conclusions. 
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Coming to the energy itself, i.e. to the integrated prod- 
uct of the interaction Vc{r) with the central energy den- 
sity, we compare 12 



FIG. 2: (Color online) Comparison of C„ n (R, r) and 
Cnn IE (R- 5 r ) where the latter is either computed from NV- 
DME or PSA-DME It-functions. Upper panels: two- 
dimensional integrands. Lower panels: ratios of C^ IE (R, r) 
over C F n (R,r) for fixed values of R. Densities are obtained 
from a self-consistent EDF calculation of 208 Pb with the SLy4 
Skynne EDF in the particle-hole part and no pairing. 



Figure |] shows 11 that both NV-DME and PSA-DME 
provide comparably good profile-reproduction of the inte- 
grand C F (R, r) within the typical range of nuclear inter- 
actions (r~2 fm). Beyond such a non locality, the quality 
of the reproduction deteriorates significantly, with that 
of PSA-DME deteriorating slightly faster. In addition, 
one sees from the lower panels of Fig. [5] that the quality 
of the reproduction decreases as one goes to the nuclear 
surface, i.e. for R > 4 fm. This could be slightly im- 
proved by taking into account the deformation of the 
local momentum distribution when designing the PSA- 
DME for the scalar part of the OBDM, which we do not 



11 Note that for semi-magic spherical nuclei used in the present 
paper, the energy densities C F n (R, r) and C^ 7 ^ IE (R, r) only de- 
pend on the magnitude of R. 



E^[nn] = AnjdIldrr 2 vc{r)C F n {R,r), (44) 
Eg ME [nn] = Air J dHdrr 2 v c {r)C™ Er R,r).{^) 



12 We do not analyze individual couplings of the Skyrme-like EDF 
produced through the DME (Eq. |20| in the present paper, but 
rather test the complete Fock energy provided by each of the 
terms (i.e. central, tensor, spin-orbit) of the two-nucleon in- 
teraction. We postpone to a forthcoming publication |22| the 
analysis of the EDF couplings computed from realistic two- and 
three-nucleon chiral interactions using appendix [B| 
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Gaussian + NV-DME 
Gaussian + PSA-DME 
Yukawa + NV-DME 
Yukawa + PSA-DME 




FIG. 3: (Color online) Percentage error of Eg ME [nn] com- 
pared to Ec[nn], where the former is either computed from 
NV-DME or PSA-DME Il-functions. Densities are obtained 
from self-consistent EDF calculations using the SLy4 Skyrme 
EDF in the particle-hole channel and no pairing. 



Figure [3] shows the relative error obtained from the two 
DME variants compared to the exact Fock contribution 
for both the Gaussian and the renormalized- Yukawa ra- 
dial form factors and for three semi-magic isotopic chains. 

Let us start with Fig. [4] that shows that the depen- 
dence of the accuracy on the range of the (Gaussian) 
interaction used is significant, i.e. about a factor of two 
between a = 1.0 fm and a = 1.5 fm. As can be expected 
from the two-dimensional density profiles in Fig. [2j the 
accuracy decreases as the range of interaction increases, 
which holds for all available DME techniques [TBI I23 - 
|2"5] . This stresses that the local quasi-separability of the 
OBDM with respect to r and R underlining the DME, 
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FIG. 4: (Color online) The same as Figure [3] but for two 
different values of the range of the Gaussian interaction. 



which is exact in INM, deteriorates with increasing non- 
locality r in finite nuclei. As long as the hypothesis of 
quasi-separability is well realized within the range of the 
interaction, the DME can be quantitatively successful. 

On average, the error obtained with PSA-DME and 
NV-DME are similar as can be seen in Fig. [3J i.e. about 
6 — 8% for the three isotopic chains and for both for the 
Gaussian and the renormalized- Yukawa interactions. In 
a future publication, we demonstrate that one can obtain 
a better accuracy (1-2% error) by using a parameterized 
and empirically optimized phase-space distribution that 
takes the diffuseness of the Fermi surface into consider- 
ation. A similar improvement over that of Ref. |16j is 
reported in Refs. [27l|28]. 



C. Fock contribution from Vt 

We now turn to the Fock contribution coming from the 
tensor part of the two-nucleon interaction. As shown by 



Eq. (12 1, such a contribution involves bilinear products 



of non-local spin densities. As a matter of fact, two terms 
with different analytical structures emerge such that the 
exchange tensor energy-density reads 13 



T n F „(R,r) 



71, (R,r) 



= T r f„ 4 (R,r) + T^ 2 (R,r) , (46) 
^- J de r s n (ri,r 2 ) ■ s n (r 2 ,ri) ,(47) 

-J— s n ^(r 1 ,r 2 ) 
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x s n i/ (r 2 ,ri) , 



where T„ n i(R, r) also appear in the central contribution 
to the Fock energy. The two DME counterparts, which 



13 We recall that the weights of the two terms have been omitted 
in agreement with Eq. 1121 
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eventually depend on which variants of the DME is being 
adopted, read 
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INM-DME 
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[TlUkpr)] 2 J2 (-W(R) Jn, M u(R) 
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and reduce for spherical systems to 
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DME 
nn,2 



(R,r) 



Q. 



n^r)] 2 J„(R)- J„(R),(49) 

(50) 



One recovers a pattern which is seen when deriving the 
empirical Skyrme EDF from an auxiliary Skyrme effec- 
tive interaction. That is, the central part of the in- 
teraction only produces the so-called symmetric bilin- 
ear tensor terms proportional to J n .^(R) J niAt „(R) while 
(R., r) that contains asymmetric bilinear tensor 
terms proportional to J n)jU ,„(R) J„ )1/jU ,(R) solely comes 
from the tensor interaction |36j . This can be easily 
traced back to the spin-space coupling that character- 
izes the tensor operator. Since the numerical tests are 
presently carried out for spherical systems, we are only 
concerned with (R, r) and T^(R,r). For spin- 
unsaturated nuclei, T 7 f n 1 (R, r) is highly localized around 
the nuclear surface as seen in Fig. [| for 208 Pb. The 
same figure shows the progressive and significant im- 
provement that the PSA approach brings to the DME of 
the vector part of the OBDM. Within the typical range 
of nuclear-interactions, NV-DME falls off much faster 
than PSA-DME. Less importantly, NV-DME also intro- 
duces artificial and pronounced structures in a region 
that corresponds to the tail of the interaction. Both of 
these drawbacks are rectified progressively by PSA-DME. 
While most of the improvement is already brought by the 
spherical PSA (P 2 (R) = 0), an even better accuracy is 
obtained by incorporating the quadrupolar deformation 
P2(R) of the local momentum Fermi distribution. The 
overestimation of T E n i (R, r) at very small r seen for all 
DMEs in the lower panels of Fig. [5] corresponds to a re- 
gion where the integrand is small and where its weight is 
further reduced in the integrated energy by the r 2 phase- 
space factor. 

Coming to the energy itself, i.e. to the integrated prod- 
uct of the interaction vt (r) with the tensor energy den- 
sity, we compare 

E^[nn] = ATrJdR,drr 2 VT(r)T^ n (R,r), (51) 
E$ ME [nn] = 47ry"dRdrr 2 VT (r)r^ Mis (R,r).(52) 

which for spherical nuclei reduce to the contribution 
from T E nl and T^[ E . Figure [6] shows the relative er- 
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FIG. 5: (Color online) Comparison of T„ nl (R,r) and 
T^ B (R,r) where the latter is computed from NV-DME, 
PSA-DME or from PSA-DME with P 2 n (R) = which we 
denote as INM-DME. Upper panels: two-dimensional inte- 
grands. Lower panels: ratios of T^ l A { E (R, r) over T^ n l (R,r) 
for fixed values of 7?. Densities are obtained from a converged 
self-consistent calculation of 208 Pb with the SLy4 Skyrme 
EDF in the particle-hole channel and no pairing. 
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compensating for its initial overestimation. In contrast, 
PSA-DME stays close to the exact value for a much larger 
range of r values. 

There exist short sequences of isotopes for which the 
percentage error shows a considerable increase. The fact 
that both DMEs display such a feature suggests that 
the problem is independent of the specific form of the 
Ilf function used. To identify the source of the prob- 
lem, Fig. FJ shows T^ nl (R,r) for three nuclei displaying 
a sudden loss of accuracy. One notices that T J ^ tl (R, r) 
extends over larger intervals in R and r than for 208 Pb 
(see Fig. [5]). This corresponds to the fact that the selected 
nuclei are nearly spin-saturated and generates very small 
E^[nn] in absolute value, as seen from the lower pan- 
els of Fig. [7j As a result, the relative inaccuracy of any 
DME becomes large and the percentage error increases 
suddenly. Of course, the resulting error in the total EDF 
remains very small as the corresponding tensor contribu- 
tion is anyway negligible, i.e. the local spin-orbit density 
J 9 (R) is close to zero in nearly spin-saturated nuclei. 
Eventually, those sudden losses of relative accuracy are 
not as worrying as Fig. [6] initially suggests. 

In conclusion, the use of PSA techniques has allowed 
us to bring the DME applicable to the bilinear product 
of non-local spin densities on the same level of accuracy 
as for terms depending on the scalar part of the OBDM. 
One could certainly work even harder to bring the overall 
DME accuracy below 1%. This could be achieved (i) by 
allowing free parameters in the II— functions to be op- 
timized on a set of reference calculations 14 and/or (ii) 
by going to higher orders in the DME, consistently for 
both the scalar and the vector parts of the OBDM. This 
should however be done within the frame of the general- 
ized Skyrme EDF proposed in Ref. [B]. 



FIG. 6: (Color online) Percentage error of E^ AIE [nn] com- 
pared to E E [nn] where the former is either computed from 
NV-DME or from PSA-DME. Densities are obtained from 
self-consistent EDF calculations using the SLy4 Skyrme EDF 
in the particle-hole channel and no pairing. Notice the differ- 
ent vertical scale compared to Fig. [3] 



ror of NV-DME and PSA-DME compared to the ex- 
act Fock contribution, for both the Gaussian and the 
renormalized- Yukawa radial form factors and for three 
semi-magic isotopic chains. For both types of interac- 
tion, the percentage error of NV-DME easily reaches 
40%. This is in contrast to PSA-DME whose percent- 
age error is typically within ±10% for most parts of the 
three isotopic chains. This can be traced to the fact that, 
while both NV-DME and PSA-DME overestimate the 
reference quantity for small r (typically less than 1 fm), 
NV-DME decreases much faster with r, thereby over- 



D. Fock contribution from Vls 

1. Basic analysis 

We now turn to the spin-orbit contribution to the Fock 
energy. As shown in Eq. ( 11 ), and unlike for central and 
tensor forces, such a contribution involves both the scalar 
and the vector parts of the OBDM. In this case, we first 
compare the spin-orbit energy density 

^l(R,r) = ~ Jde r s n (r 1 ,r 2 ) ■ r x V 2 p n (r 2 , riQ5?) 



14 As will be shown in a future publication, parameterizing II| can- 
not remove the sudden loss of relative accuracy discussed above 
for spin-saturated nuclei. As already stated, this is not a prob- 
lem in the end as the corresponding contribution to the energy 
is negligible anyway. 
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FIG. 7: (Color online) A few representative nuclei with dif- 
fuse T^ ljl (R, r) together with absolute E F [nn\ for the corre- 
sponding isotopic chains. Densities are obtained from a self- 
consistent EDF calculation using the SLy4 Skyrme functional 
in the particle-hole part and no pairing. 



to its DME counterpart 
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LS™»(K,r) = -ni(k n F r)r 



xJ„(R).V R IIg(fc£rK(R) (54) 



Note that terms containing more than two gradients have 



been truncated in LS 



DME 
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which eventually depends on which variants of the DME 
is being adopted 15 and that reduces for spherical systems 



15 The numerical tests shown in the present section actually use 
INM-DME rather than PSA-DME, i.e. k q F is employed rather 
than k q F in Hf. We still label the results as PSA-DME as no 
significant difference is seen compared to INM-DME. 



FIG. 8: (Color online) Comparison of LS™(R, r) and 
LS„^ IE (R.,r) where the latter is computed from either 
NV-DME or PSA-DME. Upper panels: two-dimensional 
integrands. Lower panels: ratios of LS^ IE (R, r) over 
LS F n (R, r) for fixed values of R. Densities are obtained from 
a converged self-consistent calculation of 208 Pb with the SLy4 
Skyrme EDF in the particle-hole channel and no pairing. 

Figure [8] shows that PSA-DME significantly over- 
estimates (in absolute values) the maximum peak of 
LS nn (R,r) at the nuclear surface. In addition, oscil- 
lations at larger r, i.e. in the tail of the two-nucleon 
interaction, are not captured by PSA-DME. In contrast, 
NV-DME reproduces relatively well the density profile 
LS E n (R, r), in particular as for the main peak at the nu- 
clear surface. This suggests that the significant improve- 
ment for PSA-DME over NV-DME as to reproducing the 
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tensor energy density does not transpose to the spin-orbit 
energy density. The previous assertions are supported 
by tests carried over several isotonic and isotopic chains. 
Looking for possible improvements, we tested that in- 
cluding truncated higher-order terms associated with the 
action of Vr on (l/4AjO„ — r„ + 3/5fc^ 2 p„), when going 
from Eq.[53jto[54] does not improve the accuracy of PSA- 
DME. 

Coming to the energy itself, i.e. to the integrated prod- 
uct of the interaction u_LS , ( r ) with the spin-orbit energy 
density, we compare 

El s [nn) = 4tt J dlidrr 2 v LS (r) LS* n (R,r) , (55) 

] = 4n J dRdrr 2 v L s(r)r 2 LS™ E CR,iJ>6) 



E 



DME 
LS 



Figure [9] shows the percentage error obtained for three 
isotopic chains. In agreement with the analysis done 
for the spin-orbit energy density, the percentage error 
of PSA-DME is impractically large and negative, in the 
range of -15% to -50% for the two schematic interactions 
used. In contrast, NV-DME provides a much better ac- 
curacy with percentage errors within ± 10% for most 
studied isotopes. Last but not least, one notes that the 
spikes in the percentage errors already discussed in sec- 
tion |IV C| arise for the same isotopes and relate to the 
vanishing non-local spin density in near spin-saturated 
nuclei. 



2. Further investigation of the spin-orbit exchange 

The results of the previous section show that NV-DME 
is better suited than PSA-DME to reproduce the spin- 
orbit contribution to the Fock energy. This can be con- 
founding in light of the better accuracy obtained using 
PSA-DME to reproduce the tensor contribution to the 
Fock energy. We can infer from Fig. [5] that NV-DME 
underestimates the main peak of the nonlocal spin den- 
sity while the latter is well captured by PSA-DME. It is 
thus puzzling to find the opposite for the Fock spin-orbit 
energy density. In the following we employ a toy model 
of the OBDM of finite nuclei to show that this is due to 
a fortuitous cancelation of errors. 

Having already a handle on the non-local spin den- 
sity s g (ri,r2), we focus on the term it multiplies in the 
spin-orbit energy density, i.e. r x V2p g (ri,r2), which we 
first approximate by r x VR/5 9 (ri, T2) thanks to the weak 
dependence of the non-local matter density on the ori- 
entation of r [31]. Hence, and focusing arbitrarily on 
neutrons, we want to compare the two quantities 



G E = V R p„(R,r) , 
Gdme = V R (ng(fc£r)p„(R) 



(57) 
(58) 



where the latter is independent of whether NV-DME or 
PSA-DME is used. To do so, we employ a toy model in 
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FIG. 9: (Color online) Percentage error of E£g IE [nn] com- 
pared to E^s[nn] where the latter is either computed from 
NV-DME or from PSA-DME. Densities are obtained from 
self-consistent EDF calculations using the SLy4 Skyrme EDF 
in the particle-hole channel and no pairing. Notice the differ- 
ent vertical scale compared to Figs. [3] and [6] 



which the nonlocal and local matter densities are built 
from a three-dimensional harmonic oscillator model with 
smeared occupancy 37J. The corresponding analytical 
expressions, as given in Ref. [37], read as 
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r r 

p n (K+ 2' R_ 2^ = 6XP 



-l/4a 



1 - t 



Pn(R) 



2cC 
Z372 



(1 - i 2 )- 3 / 2 exp 



-c^i? 2 



1 - t 



1 + t 



where a 2 = mw/K, and from J p n (R) dR = N, we have t = 1 — (2/iV) 1 / 3 . From Eqs. 
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r r 

V R p„(R+ -,R- -) = exp 



-l/4a 



■„ 2 1 + * 
1 -t 



V fl p„(R) 
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The corresponding PSA-DME reads 



Pn (R+-,R- 2 )^3— — 



1 + ^ - 
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(59) 
(60) 

and 60 one easily obtains 
(61) 
(62) 

p n (R) , (63) 



1-t 

T+t 



such that, given the definition of k F (R), one can easily 
obtain 
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U p (k F r)p(K) 



:jo(A£r)V R p„(R) (64) 



and show that 



G ra ti (R, r) = 
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DME 



(R,r) 



G E (R,r) 



j (k F r)exp 



l/4a 2 r 



2 2 1+* 
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In order to study G ra uo quantitatively, we fix the inverse 
oscillator length, a, using the Blomqvist and Molinari 
formula [3S], i.e. 1/a 2 = (0.90 A 1 / 3 + 0.70). In sub- 
sequent discussions, we take reasonable combinations of 
A and N although we show that the conclusions of the 
present section are independent of the actual value of A. 

Before analyzing the behavior of G ra u (R, r), it is 
worth noticing that the toy nonlocal matter density is 
exactly separable in relative and center-of-mass coordi- 
nates. Such a separability being one inherent, usually 
only approximate, aspect of the DME, we expect the lat- 
ter to work well in the present case [37]. Computing the 
same ratio as in G ra u (R, r) without the gradient oper- 
ators, we do indeed obtain the good performance of the 
DME as is visible in Fig. [TUJ Note in particular that 
the ratio is independent of the value of R. Such a result 
proves that the toy model provides a situation compa- 
rable to the one studied in Sec. |IVB| i.e. the DME of 
the scalar part of the density matrix performs well. Such 
a performance sets the stage in view of qualifying the 
results obtained below for G rat i (R, r). 

In order to identify the short distance behavior of 
G ra tj (R, r), we perform a Taylor series expansion in r 



GratioiR, r) « 1 



un2 



a 2 (l + t) 
6 ' 4(1 - t) 



(65) 



Looking close to the surface of the nucleus, one can 
neglect k F 2 /6 in comparison with the second term of 
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FIG. 10: (Color online) Ratio of the DME (Eq.(63l) over the 



exact (Eq.(|59[)) expressions of the toy nonlocal matter density. 



obtains 



Eq. (j65J) . Defining G error (R,r) = G rafM (R, r) - 1, one 

(66) 



a 2 (l + i) 



4(1 -t) 



the nuc 
of Eq 



Eq. (661 is valid around the nuclear surface. Inside 
eus, one cannot neglect the first term (k F 2 /6) 
( 65 ) . This is irrelevant as the spin-orbit energy 
density is concentrated around the nuclear surface. Fig- 
ure [TT] bears our expectation i.e. overestimation of Ge 
by Gdme around the nuclear surface for a wide range of 
R, A and N values. It can also be seen that there is a 
gradual and systematic shift from slight underestimation 
to overestimation as one moves from inside the nucleus 
to the nuclear surface. 

Keeping the results shown in Fig. [10] as a reference, 
we conclude that the application of the gradient oper- 
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FIG. 11: (Color online) G ra ti (R,r) as a function of r for a 
selected set of (R, A, N). 



ator on the scalar part of the density matrix deterio- 
rates the quality of the DME that overestimates the ex- 
act results, in particular as one goes to the surface of 
the nucleus where the exchange spin-orbit energy density 
is maximum. Combined with the good approximation 
of the vector part of the density matrix, such a semi- 
quantitative analysis explains the overall overestimation 
(in absolute value) of the exchange spin-orbit energy pro- 
vided by PSA-DME (see Fig. [9). Contrarily, the underes- 
timation of the vector part of the density matrix by NV- 
DME provides a fortuitous, but rather accurate, cance- 
lation of errors such that the nonlocal spin-orbit energy 
density is much better reproduced overall (see Fig. [9]). 
Even though we can be satisfied with such a situation in 
the short term future and advocate the use of the NV- 
DME variant for the spin-orbit contribution to the Fock 
energy, it would be more satisfying on the long run to 
design a suitable DME for the gradient of the scalar part 
of the density matrix that can be combined with the im- 
proved PSA-DME for the vector part. 



V. CONCLUSIONS AND OUTLOOK 

The present paper is part of a long-term project to 
build non-empirical nuclear energy density functionals 
from realistic two- and three-nucleon interactions using 
many-body perturbation theory [SHU] ■ The density ma- 
trix expansion is an important component of this effort, 
as it can be used to construct numerically-tractable ap- 
proximations to the non-local Hartree-Fock energy. In 
the first part of this paper, we assessed the accuracy of 
the DME at reproducing central, tensor, and spin-orbit 
contributions to the non-local Fock energy. Our central 
finding is that the conventional DME of Negele and Vau- 
therin performs very poorly in describing the spin-vector 
part of the density matrix, while the scalar part is de- 
scribed reasonably well. In order to address this defi- 
ciency, we have reformulated the density matrix expan- 



(i) It allows one to design expansions of both the scalar 
and the vector parts of the OBDM on an equal 
footing. This constitutes a significant improvement 
over the formulation of Negele and Vautherin who, 
as they acknowledged in their seminal paper, were 
not able to provide a satisfactory expansion of the 
vector part of the density matrix. Considering that 
the vector part of the density matrix is non-zero in 
spin-unsaturated nuclei, i.e. in the large major- 
ity of nuclei, such an improvement is mandatory 
in view of constraining a universal energy density 
functional. 

(ii) By construction, the PSA formulation allows one to 
incorporate information about the local momentum 
distribution of the Fermi system of interest. For the 
scalar part of the OBDM, one recovers the satisfac- 
tory expansion of Negele and Vautherin by averag- 
ing over the phase space of the locally- equivalent 
infinite nuclear matter system. For the vector part 
of the OBDM, one can go beyond this by taking 
into account the anisotropy that characterizes the 
local-momentum distribution at the spatial surface 
of finite Fermi systems. In contrast to the scalar 
part of the density matrix for which it has little 
impact, incorporating the deformation of the lo- 
cal momentum distribution in the expansion of its 
vector part is crucial since the latter peaks at the 
nuclear surface where such an anisotropy is maxi- 
mum. 

In the second part of the paper, we gauged the accu- 
racy of the new PSA-DME and the original NV-DME 
over a large set of semi-magic nuclei using two non-self 
consistent measures, i.e., the Fock energy density pro- 
file and the Fock energy itself. The different analytical 
structures of the central, tensor and spin-orbit contribu- 
tions led us to perform separate tests for each type of 
contribution. The main conclusions were: 

(a) A few percent accuracy is reached for the central 
force contribution to the Fock energy that depends 
on the scalar part of the density matrix. The level 
of accuracy is insensitive to the particular variant 
of density matrix expansion. 

(b) For Fock energy contributions from the central and 
tensor forces that depend on the vector part of the 
density matrix, the original expansion of Negele 
and Vautherin leads to about 50% errors. The new 
expansion based on phase-space averaging tech- 
niques reduces errors to the few percent level, which 
is the same level of accuracy as for terms involving 
the scalar part of the density matrix only. 

(c) The spin-orbit exchange is somewhat trickier as it 
combines the vector part of the density matrix with 
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the gradient of its scalar part. Surprisingly, the ex- 
pansion of Negele and Vautherin is shown to work 
much better than the new one proposed here. Using 
a semi-realistic toy model, we demonstrated that 
this is due to a fortuitous cancelation of errors be- 
tween the underestimation of the vector part of the 
density matrix and the overestimation of the gra- 
dient of its scalar part. Even though one can be 
satisfied in the short term with using the NV-DME 
variant for the spin-orbit contribution to the Fock 
energy, the present analysis calls for the design of 
a suitable expansion of the gradient of the scalar 
part of the density matrix that can be combined 
with the improved expansion proposed here for the 
vector part. 

Optimizing the density matrix expansion for the cen- 
tral, tensor and spin-orbit contributions to the Fock en- 
ergy as explained above, one reaches an overall error level 
of a few-percent over a representative set of semi-magic 
nuclei. With such an accuracy at hand, one can envision 
using the corresponding generalized Skyrme-like energy 
functional as a microscopically-constrained starting point 
around which future refined phenomcnological parame- 
terizations can be built. Indeed, the goal of a forthcom- 
ing publication 22J is to explicitly compute and analyze 
all the density-dependent couplings entering the general- 
ized Skyrme-like energy density functional starting from 
realistic two- and three-nucleon Chiral-EFT potentials at 
N 2 LO [Ml EH]- Of particular interest will be the analy- 
sis of (i) the importance of building explicit pion physics 
into the energy functionals, (ii) the density dependence 
of spin-orbit and tensor couplings in view of their anal- 
ysis in recent phenomenological studies [Ml [HJ 132] 
and (iii) the role of three-nucleon forces in these aspects, 
as well as their effects on the evolution of nuclear shells 
with isospin. Still, the EDF obtained in this approach 
will only contain the Hartree-Fock physics such that fur- 
ther correlations must be added to produce any reason- 
able description of nuclei. In the short term, such an 
addition can be done empirically by adding the DME 
couplings to empirical Skyrme functionals and perform- 
ing a refit of the Skyrme constants to data. While this 
is a purely empirical procedure, it is motivated by the 
well-known observation that a Brueckner G-matrix dif- 
fers from the vacuum NN interaction only at short dis- 
tances. Therefore, one can interpret the refit to data as 
approximating the short-distance part of the G-matrix 
with a zero-range expansion thru second order in gradi- 
ents. Eventually though, it is the goal of a future work to 
design a generalized DME that is suited to higher orders 
in perturbation theory. 

In addition to using the results of the present and forth- 
coming papers as building blocks for a microscopically- 
constrained Skyrme phenomenology, additional work is 
needed to validate the density matrix expansion method 
and to gauge its accuracy. Given the outcome of our 
analysis, several paths can be followed: 



(i) The conclusions reached in the present work must 
be further validated through self-consistent tests, 
i.e. binding energies, radii and single-particle en- 
ergies must be benchmarked against self-consistent 
Hartree-Fock calculations. The question of whether 
the Hartree term must be treated exactly is to be 
addressed quantitatively in such a context. 

(ii) An even better accuracy could be reached for the 
central and tensor contributions to the Fock energy 
by going consistently to higher orders in derivatives 
in the expansion of both the scalar and the vector 
parts of the density matrix. This should be done 
within the frame of the extended Skyrme energy 
density functional proposed in Ref. [§]. 

(iii) As already stated, the present analysis of the spin- 
orbit contribution calls for a suitable expansion of 
the gradient of the scalar part of the one-body den- 
sity matrix. 
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Appendix A: Local densities 

Non-zero local densities can be formed by taking 
derivatives of the OBDM up to second order. In the 
basis from which |$) is built, they read 



r q (r) 
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^V^J(rg) • Vtpiirq) p%, 
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■ V M y>t(rg) fi(rq))pl 
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4, (A6) 



F q Ar) = \ J2 



V M v?i(rg) 



V-a^(rg) p?. . (A7) 



and denote the matter density, the kinetic density, the 
spin density, the current density, the spin-current pseu- 
dotensor density, the spin kinetic density and the tensor 
kinetic density. In the above formulae, ipi(rq) denotes a 
spin 1/2 spinor. Among those local densities, the time- 
odd ones |24j vanish in time-reversal invariant systems, 
viz, 



s 9 ( r ) 
T,(r) 



jg(r) = 0, 
F,(r)=0. 



(A8) 



Appendix B: Skyrme-like couplings 



from the application of the DME to the Fock contribution 
to the ground-state energy (Eqs. (fl0|), (11) and (12) 



with the proper coefficients restored). The central, spin- 
orbit and tensor parts of the two-nucleon interaction are 



as specified in section II C These couplings are derived 
under the assumption of time-reversal invariance. The 
case where time-reversal invariance is relaxed will be the 
subject of a future publication. 

Starting from the definitions 



af T [nf /s nf /s ] = 4nJdrr 2 VT s (r)Ilf s Ilf s ,(Bl) 
af^nfllf] ee fj d rr% TS (r)U^ Ilf/ S ,(B2) 



We now provide explicit expressions of the couplings 
entering the Skyrme-like functional (Eq.(20)) that results 



the couplings take the form 
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To carry on further the computation of the couplings, one 
must choose an explicit form of the two-nucleon interac- 



tion and perform the integrals entering Eqs. Bl and B2 



As schematic interactions have been used in the present 



paper for illustrative purposes, we postpone such an in- 
tegration to the explicit computation of the couplings 
obtained from a Chiral-EFT lagrangian at N 2 LO 22]. 
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Appendix C: Local anisotropy P2(r) 

The Husimi distribution is one of the many quantum 
phase-space distribution functions. It possesses the key 
property of positive definitencss 32 , 43 and is defined as 



H q (r,p) 
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ipi(Tiq)e 2r o dr 1 
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where N = 1/ (tt 3 ^^ 2 ) and ro is a chosen parameter. 
To derive Eq. (22) for the quadrupolar local anisotropy 



of the momentum Fermi surface P^ (r) we start from the 
definition 
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and make use of the relations 
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Through direct application of the above relations, one 
obtains 

fdpp 2 H q (r lP ) » (2^) 3 ^ 5 ^|V^(r (7 )| 2 p« 
+O((4r ) 2 ) , 
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which, plugged into Eq.( |C2[ ), gives 
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harmonic relations 
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turn out to be handy. In these relations, Y™ 1 refers to 
a spherical harmonic function and Pi refers to Legendre 
polynomial of order 1. Applying these relations, one ob- 
tains 
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where F(l,j) is some function of i and j. The occupa- 
tion probability of a given spherical shell p qn ^ 1 is one or 
zero, except for open-shell semi-magic nuclei where the 
so-called filling approximation provides the valence shell 
with a partial occupation. To obtain the explicit form of 
F(l,j), one can use the relation 
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Further simplifications can be performed for spherical 
systems, using single-particle wave- functions expressed in 
terms of spherical coordinates r = (r, 9, if) as 



and perform involved angular momentum coupling oper- 

ations. Alternatively, one notes that X^l^ 7 ^ 1 "?)! P%. IS 
nothing but the kinetic energy density given in Eq. ( A2 ) 
and use the corresponding expression (44] . Either way, 
one obtains 



<Pi(rq) 



l(l + l)(2j + l) 

47T 



(C8) 



through several angular momentum coupling operations. 
For that, the following Clebsch-Gordan and spherical 



Plugging these intermediate results into Eq. ( C5 ) yields 
the expression of P2C 1 *) a s 



,(r) 



E 

nlj 



2j + 1 

47T 



dr r 



(C9) 



22 



where 



r q (r) 



E 



2.? + 1 

47T 



8 V^(r) 
dr r 



(CIO) 



Appendix D: Scalar part of the OBDM 



Appendix E: Vector part of the OBDM 



We start from Eq. (24), average over the orientation 
of k and r 16 , and apply relations 



l - J de r (r-A)(r-B) = r - A ■ B , (Dl) 
= V 2 /?(r)-2r(r), (D2) 



4tt 

V 2 + V' 2 ]p(r,r') 



to obtain 

p 9 (R 



r r N 

2' R ~2' 



j (kr) p q (R) + L(kr) p q (R) (D3) 
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jo(fcr) (Ap 9 (R)-4r,(R)) 



where 



f/fcr) 2 

L(fcr) = 2fcr jx(fcr) - i— i- j (fcr) . 



(D4) 



As discussed in section III E| the effects of anisotropy 
and diffuseness are minimal for the scalar part of the 
OBDM. Therefore, we perform the PSA over the phase- 
space of the locally equivalent pure-isospin nuclear mat- 
ter 17 to obtain 



p q (R+l,R-l) 



Ji(k q F r ) 
k q F r 
r 2 ji{k q F r) 



1 V,(R) 

F' 

2 A. (VI , , 

g q (R), (D5) 



2 k q F r 

with the second-order correction density being composed 
of 



0q {R) = ~Ap g (R) - r q (R) + \kf A(fc«r)p,(R) .(D6) 



Expanding A(k F r) in Taylor series, one has 
A(k q F r)^l + 0((k q F r) 2 ). 



(D7) 



such that, by retaining the lowest order only, one recovers 
Eq. [26] with the II— functions given by Eqs. [27] and [28[ 



16 The order of the two averaging operations is dictated only by 
the requirement of simplicity. In this case, we averaged over the 
orientation of r followed by that of k. 

17 The angle integration with respect to the orientation of k is triv- 
ial as such a dependence has already been averaged out. 



We start from Eq. (34). For time- reversal invariant 



systems, the local spin density s ? (r) vanishes. Conse- 
quently, the only non-vanishing contribution relates to 
the term r ■ (Vi — V2). Using the definition for the lo- 
cal spin-current pseudotensor density given by Eq. ( A5 1 , 
one obtains 



R 



;>R 



k \ " 



r M J^uiR) -(El) 



ft 



The final step involves performing the PSA over a de- 
formed sphere that characterizes the local momentum 
distribution. Let us start from a spheroid given in mo- 
mentum space given by the equation 



k' 2 , 



a(R) 2 a(R) 2 c(R) 2 



1 . 



(E2) 



For ease of notation, we write a(R) as a and c(R) as 
c in the following. We constrain the position-dependent 
quantities a and c by requiring that the spheroid has a 
given volume and quadrupole moment, viz, 



V q ee ^ 3 = ^a 2 c, 



2(-a 2 + c 2 ) 
2 a 2 + c 2 



(E3) 
(E4) 



The II— function is obtained via the integration over the 
phase space of interest 



4w*k F 3 Jv q 



dk e 



-irk 



(E5) 



Carrying out the integration over the volume V q encom- 
passed by the spheroid given in Eq. (E2) can be done 



by using a stretched coordinate system from the trans- 
formation 



k {kxi ky,kz) y k (k x , ky, —kg). 

such that one finally obtains 



(E6) 



9)^ 



fi—X 



where 



Ul(k q F r 



. h{k q F r) 
k F r 



(E8) 
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and In spherical systems, both the pseudoscalar and the ten- 

sor parts vanish such that one obtains 



^2 P|(R) ) ■ [ yj 

Setting P| = Oj which consists of performing the PSA 
over INM phase-space, results in the same II— function 
with k q F replaced by k q F . 

For spherical systems, one can simplify the expression 
further by writing J 9j(La ,(R) as a sum of pseudoscalar, 
vector and (antisymmetric) traceless tensor parts 



•WR) = l^J^W + l E^ fc jW(R) 

+J(^(R), (ElO) sJR+l,R-lj c --^(fc^rx J 9 (R). (E14) 



where the three components read 

z 

J(°)(R) = 2 <5^J^(R), (Ell) 

(J,,I/=X 

z 

J 2( R ) = E ^-W*). (E12) 

fJ,,V=X 

4EW$(R). (E13) 
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